import std.math;
import std.stdio;
import std.conv;


    // compute the FFT of x[], assuming its length is a power of 2
    cdouble[] fft(cdouble[] x) {
        long N = x.length;

        // base case
        if (N == 1) return [ x[0] ];

        // radix 2 Cooley-Tukey FFT
        if (N % 2 != 0) { throw new Exception("N is not a power of 2: "~to!string(N)); }

        // fft of even terms
        cdouble[] even = new cdouble[N/2];
        for (long k = 0; k < N/2; k++) {
            even[k] = x[2*k];
        }
        cdouble[] q = fft(even);

        // fft of odd terms
        cdouble[] odd  = even;  // reuse the array
        for (long k = 0; k < N/2; k++) {
            odd[k] = x[2*k + 1];
        }
        cdouble[] r = fft(odd);

        // combine
        cdouble[] y = new cdouble[N];
        for (long k = 0; k < N/2; k++) {
            double kth = -2 * k * PI / N;
            cdouble wk = cos(kth)+ sin(kth)*1.0i;
            y[k]       = q[k]+(wk*(r[k]));
            y[k + N/2] = q[k]-(wk*(r[k]));
        }
        return y;
    }


    // compute the inverse FFT of x[], assuming its length is a power of 2
    cdouble[] ifft(cdouble[] x) {
        long N = x.length;
        cdouble[] y = new cdouble[N];

        // take conjugate
        for (long i = 0; i < N; i++) {
            y[i] = conj(x[i]);
        }

        // compute forward FFT
        y = fft(y);

        // take conjugate again
        for (long i = 0; i < N; i++) {
            y[i] = conj(y[i]);
        }

        // divide by N
        for (long i = 0; i < N; i++) {
            y[i] = y[i]*(1.0 / N);
        }

        return y;

    }
    
    double[] fftMagnitudes(short[] sound){
    	cdouble[] dat = new cdouble[sound.length];
    	foreach(i,v; sound){
    		dat[i] = v + 0i;
    	}
    	//writeln(sound);
    	cdouble[] freq = fft(dat);
    	double[] mag = new double[freq.length/2];
    	foreach(i, ref v; mag){
    		v = abs(freq[i]);
    	}
    	return mag;
    }
/+    
    double[] freqMaximums(short[] sound, double ratePerSecond, int span){
    	double[] mag = fftMagnitudes(sound);
    	//writeln(`magnitudes: `,mag);
    	double[] ret;
    	foreach(i; span..mag.length-span){
    		bool ismax=true;
    		foreach(j; -span..span ){
	    		if( mag[i] < mag[i+j] ){
	    			ismax=false;
	    		}
	    	}
    		if( mag[i]!=0 && ismax ){
	    	  ret ~= ratePerSecond / sound.length * i;
	    	}		
    	}
    	return ret;
    }
+/    
    void main_fft_test(){
    	cdouble[] x = new cdouble[16];
    	foreach(i,ref v; x){
//    		v = sin(2*PI*i/x.length * 2) + 0i;
    	    v = 0+0i;
    	}
    	x[3] = 5+0i;
    	cdouble[] freqx = fft(x);
    	foreach(i,v; freqx){
    		writeln(i,' ',v);
    	}
    	
    	writeln();
    	cdouble[] undo = ifft(freqx);
    	foreach(i,v; undo){
    		writeln(i,' ',v);
    	}
    	
    	short[] sound = new short[1024];
    	foreach(i,ref v; sound){
    		v = cast(short)(30000*sin(2*PI*i/sound.length * 32));
    	}
    	//writeln(freqMaximums(sound, 1024));
    }
